Modelos lineales

Regresión y factoriales (ANOVA y ANCOVA)

Introducción

  • Definición

  • Tipos

  • Supuestos

  • Flujo de trabajo

Definición

Relación lineal entre dos o más variables, dependiendo una (variable dependiente) de la(s) otra(s) (variable(s) explicativa(s))

Variable dependiente: y

Variable(s) explicativa(s): x

Modelo matemático lineal: y = 2x o y = 3 - 2x

Objetivos

  1. Descripción de la relación

  2. Inferencia o generalización de los resultados

  3. Predicciones

Tipos

Variable dependiente: siempre cuantitativa

Variables explicativas:

  • Cuantitativas: modelos de regresión (simple o múltiple)

  • Cualitativas o categóricas: ANOVA

  • Cuantitativas y categóricas: ANCOVA

Supuestos del modelo

  1. Independencia

  2. Linealidad

  3. Normalidad

  4. Homocedasticidad

Flujo de trabajo

Pasos a seguir análisis de datos con modelos lineales (Cayuela y de la Cruz, 2022)

Regresión simple y regresión múltiple

Estructura

Ecuación regresión:

yi = β0 + β1 x1 + Ɛi

Ɛi ~ N(0,σ2)

β0: ordenada en el origen

β1: pendiente

Ɛi: residuos

Ajuste modelo lineal

Función lm(y ~ x)

Ejemplo

Pregunta: Masa corporal especies pingüinos y cómo esta es influenciada por la longitud de la aleta de los pingüinos

Hipótesis

H0 (hipótesis nula): no existe relación entre la masa corporal de los pingüinos y la longitud de sus aletas

H0 : β1 = 0

H1 (hipótesis alternativa): pendiente de la recta que relaciona masa corporal con longitud aleta es distinta de cero

H0 : β1 ≠ 0

1. Leer los datos en R y explorar gráficamente la relación entre las variables

library(palmerpenguins)
data(penguins)
head(penguins)
# A tibble: 6 × 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie  Torgersen           39.1          18.7               181        3750
2 Adelie  Torgersen           39.5          17.4               186        3800
3 Adelie  Torgersen           40.3          18                 195        3250
4 Adelie  Torgersen           NA            NA                  NA          NA
5 Adelie  Torgersen           36.7          19.3               193        3450
6 Adelie  Torgersen           39.3          20.6               190        3650
# ℹ 2 more variables: sex <fct>, year <int>
# eliminamos NAs variables que vamos a utilizar
library(dplyr)
penguins_clean <- penguins %>%
  filter(!is.na(flipper_length_mm), !is.na(body_mass_g))
library(ggplot2)
ggplot(data = penguins_clean, aes(x = flipper_length_mm, y = body_mass_g)) +  
  geom_point(color = "blue", size = 2) +      
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (g)",                   
      title = "Relación entre masa corporal y longitud aleta "   
  ) +
  theme_minimal()

2. Ajuste modelo de regresión

lm_reg1 <- lm(body_mass_g ~ flipper_length_mm, data = penguins_clean)

3. Variación explicada por el modelo y resolución hipótesis

Tabla del análisis de la varianza

anova(lm_reg1)
Analysis of Variance Table

Response: body_mass_g
                   Df    Sum Sq   Mean Sq F value    Pr(>F)    
flipper_length_mm   1 166452902 166452902  1070.7 < 2.2e-16 ***
Residuals         340  52854796    155455                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Suma de cuadrados: variabilidad de la variable y

Coeficiente de determinación (R2): bondad de ajuste del modelo (0-1)

4. Interpretación del modelo

summary(lm_reg1)

Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = penguins_clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1058.80  -259.27   -26.88   247.33  1288.69 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5780.831    305.815  -18.90   <2e-16 ***
flipper_length_mm    49.686      1.518   32.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 394.3 on 340 degrees of freedom
Multiple R-squared:  0.759, Adjusted R-squared:  0.7583 
F-statistic:  1071 on 1 and 340 DF,  p-value: < 2.2e-16

y = 49.686x - 5780.831

Fiabilidad estimadores para predecir

R2 y R2 ajustado

5. Predicciones con el modelo

library(broom)
library(tidyverse)
library(dplyr)
preds <- augment(lm_reg1, type.predict = "response", se_fit = TRUE)
preds <- preds %>%
  mutate(lower = .fitted - 1.96 * .se.fit,upper = .fitted + 1.96 * .se.fit)
  
ggplot(data = preds, aes(x = flipper_length_mm, y = body_mass_g)) +  
  geom_point(color = "blue", size = 2) +     
  geom_line(aes(y = .fitted)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (g)",                   
      title = "Predicciones modelo de regresión"   
  ) +
  theme_minimal()

6. Revisión supuestos del modelo

Supuestos de linealidad, normalidad y homocedasticidad

par(mfrow=c(2,2))
plot(lm_reg1)

Supuestos 1

Supuestos 2

Supuestos 3

Supuestos 4

Normalidad:

shapiro.test(residuals(lm_reg1))

    Shapiro-Wilk normality test

data:  residuals(lm_reg1)
W = 0.99301, p-value = 0.1123

Linealidad:

library(lmtest)
resettest(lm_reg1)

    RESET test

data:  lm_reg1
RESET = 15.998, df1 = 2, df2 = 338, p-value = 2.3e-07

Homocedasticidad:

library(car)
ncvTest(lm_reg1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 2.177953, Df = 1, p = 0.14

No se cumplen supuestos del modelo: transformación variables (logarítmica, raíz cuadrada, Box-Cox, log-log, etc.)

Regresión múltiple

yi = β0 + β1x1 + … + β1xn + Ɛi

Ɛi ~N(0,σ2)

Suma de cuadrados tipo I o III

Suma de cuadrados: tipos (Cayuela y de la Cruz,) 2022

Suma de cuadrados tipo I: anova()

Suma de cuadrados tipo III: Anova(lm,type=“III”)(library(car))

Modelos factoriales: ANOVA y ANCOVA

Análisis de la varianza (ANOVA)

Hipótesis: comparación medias distintos grupso (niveles de factor)

H0: medias de los grupos son iguales

H0: µa = µb = … = µ

H1: al menos una media es distinta a las demás

H0: Ǝµj ≠ µ

Ejemplo

Pregunta: Diferencias longitud pico tres especies de pingüino

Hipótesis

H0 (hipótesis nula): no existe diferencias en las medias de la longitud de los picos para las tres especies

H0 : µadelie = µgentoo = µchinstrap = µ

H1 (hipótesis alternativa): la media de al menos una de las tres especies es distinta a las demás

H0 : Ǝµj ≠ µ

1. Explorar gráficamente la relación entre las variables

ggplot(data = penguins, aes(x = species, y = bill_length_mm, fill = species)) +  
  geom_boxplot(alpha = 0.8, color = "gray30") +
  scale_fill_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  labs(x = "Species",         
       y = "Longitud pico (mm)",                   
      title = "Diferencias longitud pico entre especies"   
  ) +
  theme_minimal()

2. Ajuste modelo factorial

penguins <- penguins %>%
  mutate(species = as.factor(species))
lm_reg2 <- lm(bill_length_mm ~ species, data = penguins)

3. Variación explicada por el modelo y resolución hipótesis

Tabla de análisis de la varianza
anova(lm_reg2)
Analysis of Variance Table

Response: bill_length_mm
           Df Sum Sq Mean Sq F value    Pr(>F)    
species     2 7194.3  3597.2   410.6 < 2.2e-16 ***
Residuals 339 2969.9     8.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Interpretación del modelo y comparaciones post-hoc

summary(lm_reg2)

Call:
lm(formula = bill_length_mm ~ species, data = penguins)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.9338 -2.2049  0.0086  2.0662 12.0951 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       38.7914     0.2409  161.05   <2e-16 ***
speciesChinstrap  10.0424     0.4323   23.23   <2e-16 ***
speciesGentoo      8.7135     0.3595   24.24   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.96 on 339 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.7078,    Adjusted R-squared:  0.7061 
F-statistic: 410.6 on 2 and 339 DF,  p-value: < 2.2e-16

Variables indicadoras o variables dummy

yi,j = β1x1 + β2x2 + β3x3 + Ɛi

Un nivel es la coordenada en el origen (nivel de referencia)

 contrasts(penguins$species)
          Chinstrap Gentoo
Adelie            0      0
Chinstrap         1      0
Gentoo            0      1

Ecuación: yi = 38.79 + 10.04xi + 8.71x2

library(multcomp)
posthoc_species <- glht(lm_reg2, linfct=mcp(species="Tukey"))
summary(posthoc_species)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = bill_length_mm ~ species, data = penguins)

Linear Hypotheses:
                        Estimate Std. Error t value Pr(>|t|)    
Chinstrap - Adelie == 0  10.0424     0.4323  23.232  < 0.001 ***
Gentoo - Adelie == 0      8.7135     0.3595  24.237  < 0.001 ***
Gentoo - Chinstrap == 0  -1.3289     0.4473  -2.971  0.00895 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

5. Representación gráfica diferencias significativas

ggplot(data = penguins, aes(x = species, y = bill_length_mm, fill = species, color = species)) +  
  geom_boxplot(alpha = 0.8, color = "black") +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2, shape=21, color="black") +
  scale_fill_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  scale_color_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
 labs(x = "Species",         
       y = "Longitud pico (mm)",                   
      title = "Diferencias longitud pico entre especies"   
  ) +
  theme_minimal() +
  geom_text(data = penguins[2,], aes(y = 65, label = "a"), color="black", 
            size=5, nudge_x = 0) +
  geom_text(data = penguins[2,], aes(y = 65, label = "b"), color="black", 
            size=5, nudge_x = 1) +
  geom_text(data = penguins[2,], aes(y = 65, label = "c"), color="black", 
            size=5, nudge_x = 2)

6. Revisión supuestos del modelo

par(mfrow=c(2,2))
plot(lm_reg2)

Cambiar nivel de referencia del factor

levels(penguins$species)
[1] "Adelie"    "Chinstrap" "Gentoo"   
penguins$species <- relevel(penguins$species, ref = "Gentoo")
levels(penguins$species)
[1] "Gentoo"    "Adelie"    "Chinstrap"

Análisis de la covarianza (ANCOVA)

Hipótesis: homogeneidad de pendientes

H0: las pendientes son iguales

H1: al menos una pendiente es distinta a las demás

Ejemplo

Pregunta: La masa corporal de los pingüinos se ve influenciada por la longitud de su pico, teniendo en cuenta el sexo del individuo

Hipótesis

H0: no hay relación entre la masa corporal y la longitud del pico

H0: descontando el efecto de la longitud del pico (si existe), no hay relación entre la masa corporal y el sexo del pingüino

Homogeneidad de pendientes

Asunción adicional a las de ANOVA y regresión: relación entre la variable respuesta y la covariable (pendiente de la regresión) es la misma en todos los grupos (homogeneidad de pendientes)

ggplot(subset(penguins, !is.na(sex)), aes(x = bill_length_mm, y = body_mass_g)) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ sex) +
  labs(
    title = "Relación longitud pico y masa corporal",
    x = "Longitud pico (mm)",
    y = "Masa corporal (g)"
  ) +
  theme_minimal()

Ajuste modelo

lm_3 <- lm(body_mass_g ~ bill_length_mm*sex, data = penguins)

Comprobación de hipótesis

anova(lm_3)
Analysis of Variance Table

Response: body_mass_g
                    Df    Sum Sq  Mean Sq  F value    Pr(>F)    
bill_length_mm       1  74792533 74792533 191.8812 < 2.2e-16 ***
sex                  1  12051930 12051930  30.9194 5.572e-08 ***
bill_length_mm:sex   1    175714   175714   0.4508    0.5024    
Residuals          329 128239489   389786                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lm_4 <- lm(body_mass_g ~ bill_length_mm+sex, data = penguins)
anova(lm_4)
Analysis of Variance Table

Response: body_mass_g
                Df    Sum Sq  Mean Sq F value    Pr(>F)    
bill_length_mm   1  74792533 74792533 192.201 < 2.2e-16 ***
sex              1  12051930 12051930  30.971 5.427e-08 ***
Residuals      330 128415203   389137                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(subset(penguins, !is.na(sex)), aes(x = bill_length_mm, y = body_mass_g, color = sex)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("#F8BBD0", "#B3E5FC")) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Relación longitud pico y masa corporal",
    x = "Longitud pico (mm)",
    y = "Masa corporal (g)",
    color = "Sexo"
  ) +
  theme_minimal()
ggplot(data = subset(penguins, !is.na(sex)), aes(x = sex, y = body_mass_g, fill = sex, color = sex)) +  
  geom_boxplot(alpha = 0.8, color = "black") +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2, shape=21, color="black") +
  scale_fill_manual(values = c("#F8BBD0", "#B3E5FC")) +
  scale_color_manual(values = c("#F8BBD0", "#B3E5FC")) +
 labs(x = "Sexo",         
       y = "Masa corporal (g)",                   
      title = "Diferencias masa corporal entre sexo"   
  ) +
  theme_minimal() +
  geom_text(data = penguins[2,], aes(y = 6500, label = "a"), color="black", 
            size=5, nudge_x = 0) +
  geom_text(data = penguins[2,], aes(y = 6500, label = "b"), color="black", 
            size=5, nudge_x = 1)

Pendientes significativamente distintas

Relación masa corporal con longitud de la aleta, teniendo en cuenta la especie

ggplot(subset(penguins, !is.na(species)), aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~ species) +
  labs(
    title = "Relación longitud aleta y masa corporal",
    x = "Longitud aleta (mm)",
    y = "Masa corporal (g)"
  ) +
  theme_minimal()
datos_limpios <- penguins %>%
  filter(!is.na(flipper_length_mm))
lm_5 <- lm(body_mass_g ~ flipper_length_mm*species, data = datos_limpios)
anova(lm_5)
Analysis of Variance Table

Response: body_mass_g
                           Df    Sum Sq   Mean Sq  F value    Pr(>F)    
flipper_length_mm           1 166452902 166452902 1211.946 < 2.2e-16 ***
species                     2   5187807   2593904   18.886 1.686e-08 ***
flipper_length_mm:species   2   1519564    759782    5.532  0.004327 ** 
Residuals                 336  46147424    137344                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm_5)

Call:
lm(formula = body_mass_g ~ flipper_length_mm * species, data = datos_limpios)

Residuals:
    Min      1Q  Median      3Q     Max 
-911.18 -251.93  -31.77  197.82 1144.81 

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        -6787.281   1124.195  -6.037 4.14e-09 ***
flipper_length_mm                     54.623      5.174  10.557  < 2e-16 ***
speciesAdelie                       4251.444   1427.332   2.979  0.00311 ** 
speciesChinstrap                    3750.085   1676.687   2.237  0.02597 *  
flipper_length_mm:speciesAdelie      -21.791      6.941  -3.139  0.00184 ** 
flipper_length_mm:speciesChinstrap   -20.049      8.190  -2.448  0.01487 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 370.6 on 336 degrees of freedom
Multiple R-squared:  0.7896,    Adjusted R-squared:  0.7864 
F-statistic: 252.2 on 5 and 336 DF,  p-value: < 2.2e-16

yi = -6787.3 + 4251.4 Adelie + 3750.1 Chinstrap + (54.6 - 21.8 Adelie -20.0) flipper_length

Gráfica de predicciones

newdata <- expand.grid(
  flipper_length_mm = seq(min(datos_limpios$flipper_length_mm), max(datos_limpios$flipper_length_mm), length.out = 100),
  species = levels(datos_limpios$species)
)
predicciones_grid <- augment(lm_5, newdata = newdata)

ggplot(datos_limpios, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point(alpha = 0.5) +
  geom_line(data = predicciones_grid, aes(x= flipper_length_mm, y = .fitted, group = species, color = species), size = 1.2) +
  scale_color_manual(values = c("#F8BBD0", "#C8E6C9", "#B3E5FC")) +
  labs(
    title = "ANCOVA: predicciones por especie",
    x = "Longitud aleta (mm)",
    y = "Masa corporal (g)",
    color = "Especie"
  ) +
  theme_minimal()

Supuestos del modelo

par(mfrow=c(2,2))
plot(lm_5)